Contents

true

This document contains all data analyses presented in the main text as well as supplementary materials mentioned and extra information on the analysis for those interested (not cross-referenced).

The analysis was organized using the make-line pipeline tool provided by targets (Landau 2021). You will find raw data as well as all functions used in the pipeline in the project’s research compendium created with the rcompendium package (Casajus N. 2022).

You can take a look at the pipeline in order to have an idea of the different steps for creating the main models and figures. Each target is defined at the _targets.R file in the research compendium.

Analysis pipeline. Total runtime and storage size of each object (target) are shown.

Table of figures

Supplementary Main text
(fig-pcafullsel?) (thm-fig-1?)
(fig-clust?) (thm-fig-2?)
(fig-corcomp?) (thm-fig-3?)
(fig-data?) (thm-fig-4?)
(fig-phylord?) (thm-fig-5?)
(fig-fam?) (thm-fig-6?)
(fig-cor?) (thm-fig-7?)
(fig-psJ?) (thm-fig-8?)
(fig-psN2?)

Table of tables

Unfortunately links only redirect to tabsets open by default but not hidden ones. Please look for tables related to analyses by phylum or by compartment in their respective tabsets.

Supplementary Tables
(tbl-S1?)
(tbl-S2?)
(tbl-S3?)
(tbl-S4?)
(tbl-rich-phyl?)
(tbl-rich-trait?)
(tbl-bdtot?)
(tbl-bdann?)
(tbl-bdart?)
(tbl-bdmol?)
(tbl-bdepi?)
(tbl-bdinf?)
(tbl-bdint?)

Packages

library(dplyr)
library(tidyr)
library(vegan)
library(ggplot2)
library(targets)
library(ggrepel)
library(patchwork)
theme_set(theme_light())

Habitat Complexity

Take a look at the structure of the final complexity dataset

tar_load(comp)
knitr::kable(head(comp))
Mini Site Point Sample Branching_density Broken_density D_bin D_gray L I S Total_Density Dry_weight Sphericity DR1 DR2 DR3
MBE_1 Belle-ile Belle-ile 1 1 5.8 0.0 1.1231 2.6562 12.85333 9.143333 8.303333 181696 183.61 0.8371338 0.6460062 0.8153846 0.7113589
MBE_1 Belle-ile Belle-ile 1 2 1.0 0.0 1.0654 2.6463 16.07667 7.883333 5.393333 181696 183.61 0.6122605 0.3354758 0.7669267 0.4903587
MBE_1 Belle-ile Belle-ile 1 3 3.0 0.2 1.0722 2.6450 15.30000 8.103333 7.190000 181696 183.61 0.7470808 0.4699346 0.8873818 0.5296296
MBE_1 Belle-ile Belle-ile 1 4 4.8 0.4 1.1078 2.6492 12.91667 8.476667 6.203333 181696 183.61 0.7057078 0.4802581 0.6613704 0.6562581
MBE_1 Belle-ile Belle-ile 1 5 4.2 0.0 1.1305 2.7019 24.54667 14.523333 8.970000 181696 183.61 0.6088477 0.3654264 0.6434838 0.5916621
MBE_1 Belle-ile Belle-ile 1 6 3.0 0.4 1.1198 2.6641 16.35667 9.416667 6.783333 181696 183.61 0.6684949 0.4147137 0.7249304 0.5757082

Most variables are not normally distributed so we apply a box-cox (Box and Cox 1964) transformation and check correlations between the transformed variables

tar_read(pairscomp)
Correlogram of box-cox transformed HC variables

Correlogram of box-cox transformed HC variables

PCA of all observations

Let’s check it without removing highly correlated variables

tar_load(pcatotal)
tar_read(sptot)
Screeplot of the PCA including all observations and all variables

Screeplot of the PCA including all observations and all variables

We had up to 6 relevant axes to look at, but we know we have several highly correlated variables

tar_read(pcatot)
PCA of all observations and all variables

PCA of all observations and all variables 

tar_read(pcatotsel)
## Call: rda(X = num_sel, scale = TRUE)
## 
##               Inertia Rank
## Total               7     
## Unconstrained       7    7
## Inertia is correlations 
## 
## Eigenvalues for unconstrained axes:
##    PC1    PC2    PC3    PC4    PC5    PC6    PC7 
## 2.1984 1.6946 0.9384 0.8237 0.6952 0.4496 0.2001
tar_read(sptotsel)
Screeplot of the PCA including all observations and selected variables

Screeplot of the PCA including all observations and selected variables

Let’s check the biplots of the selected variables

tar_read(pcasel)
PCA of all observations after variable selection

PCA of all observations after variable selection

PCA of median values

tar_read(spmedsel)
Screeplot of the PCA of median values at the point level

Screeplot of the PCA of median values at the point level

Let’s check the biplots of the selected variables

tar_read(pcamed)
Figure 1 (main text)

Figure 1 (main text)

Cluster of median site values

Let’s check if we find similar groupings using a ward’s cluster

tar_read(clustcomp)
Ward's cluster of median complexity values at the site level

Ward’s cluster of median complexity values at the site level

Correlation between the centroids of the total PCA and the scores of the reduced PCA

tar_load(corcomp)
plot(corcomp[[1]])
plot(corcomp[[2]])
tar_read(cortestcomp1)
## 
##  Pearson's product-moment correlation
## 
## data:  med_sel2$PC1_c and med_sel2$PC1_score
## t = -12.623, df = 28, p-value = 4.471e-13
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9626648 -0.8416107
## sample estimates:
##        cor 
## -0.9222483
tar_read(cortestcomp2)
## 
##  Pearson's product-moment correlation
## 
## data:  med_sel2$PC2_c and med_sel2$PC2_score
## t = -14.294, df = 28, p-value = 2.165e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9702598 -0.8722282
## sample estimates:
##        cor 
## -0.9378045
Relation between Full PCA centroids and Reduced PCA scoresRelation between Full PCA centroids and Reduced PCA scores

Relation between Full PCA centroids and Reduced PCA scores

Macrofauna data availabity

tar_load(faunadata)
knitr::kable(head(faunadata))
Point Site Replicate Method Habitat Year Season Date Species Abundance
Belle-ile 1 Belle-ile 1 Grab Maerl 2007 Spring 2007-03-11 00:00:00 Abludomelita gladiosa 1
Belle-ile 1 Belle-ile 1 Grab Maerl 2007 Spring 2007-03-11 00:00:00 Abra alba 1
Belle-ile 1 Belle-ile 1 Grab Maerl 2007 Spring 2007-03-11 00:00:00 Ampharete spp. 1
Belle-ile 1 Belle-ile 1 Grab Maerl 2007 Spring 2007-03-11 00:00:00 Amphipholis squamata 2
Belle-ile 1 Belle-ile 1 Grab Maerl 2007 Spring 2007-03-11 00:00:00 Anapagurus hyndmanni 1
Belle-ile 1 Belle-ile 1 Grab Maerl 2007 Spring 2007-03-11 00:00:00 Aonides oxycephala 3 
samplespoint <- faunadata %>% 
  mutate(Date=as.factor(format.Date(faunadata$Date,"%Y"))) %>% 
  distinct(Site, Year, Point, Replicate) %>%
  group_by(Site, Year, Point) %>%
  mutate(Year = as.factor(Year)) %>% 
  summarise(count = n()) %>%
  ungroup() %>%
  mutate(State = factor(case_when(count==3 ~ "Good", count < 3 ~ "Sample(s) missing", count > 3 ~ "Too many samples")))

ggplot(samplespoint, aes(x=Year,y=Point,fill=State)) +
  geom_tile(alpha=0.5) +
  scale_fill_manual(values=c("palegreen3","lightblue", "indianred")) +
  geom_text(aes(label=count)) +
  theme(text=element_text(size =15), axis.text.x = element_text(angle = 90, size = 15), legend.position = "bottom")
Number of grab samples taken every year at each point. Empty cells show points where no samples were collected at a given date due to field constraints.

Number of grab samples taken every year at each point. Empty cells show points where no samples were collected at a given date due to field constraints.

Local macrofaunal diversity

Total richness

# Total richness
tar_load(faunaclass)
total_rich <- nrow(faunaclass %>% 
  distinct(., Species)) #725 species

# Total richness by sediment position

epi_rich <- nrow(faunaclass %>% 
  filter(Position == "Epifauna") %>% 
  distinct(., Species)) #341 epifaunal species


inf_rich <- nrow(faunaclass %>% 
  filter(Position == "Infauna") %>% 
  distinct(., Species)) #266  infaunal species


int_rich <- nrow(faunaclass %>% 
  filter(Position == "Interstice") %>% 
  distinct(., Species)) #118 interstitial species

There are 725 species in total, 341 being epifaunal, 266 being infaunal, and 118 interstitial.

Relative abundance

tar_read(relab_phy)
tar_read(relab_ord)
Average Relative Abundance of Main Phyla and OrdersAverage Relative Abundance of Main Phyla and Orders

Average Relative Abundance of Main Phyla and Orders

p.s.: it is normal to have NAs for some orders as higher classification for some of the species is ambiguous or unknown.

tar_read(relab_arth)
tar_read(relab_ann) 
tar_read(relab_mol)
Families Average Relative Abundance for Each Main PhylaFamilies Average Relative Abundance for Each Main PhylaFamilies Average Relative Abundance for Each Main Phyla

Families Average Relative Abundance for Each Main Phyla

Effects of HC in local diversity

tar_read(correl)
Correlogram of HC and physical environmental variables

Correlogram of HC and physical environmental variables

Density of macrofauna

Total macrofauna

tar_read(densselR2)
##      variables order          R2     R2Cum  AdjR2Cum         F pvalue
## 1 Current_mean     6 0.176946460 0.1769465 0.1745677 74.385774  0.001
## 2       Gravel     2 0.040422852 0.2173693 0.2128323 17.819240  0.001
## 3    PC1_score     9 0.054967407 0.2723367 0.2659908 25.985629  0.001
## 4       T_mean     4 0.021725307 0.2940620 0.2858295 10.555857  0.003
## 5    Fetch_max     8 0.014197289 0.3082593 0.2981461  7.019209  0.007
## 6           OM     3 0.029000418 0.3372597 0.3255986 14.921596  0.001
## 7        Depth     7 0.008082616 0.3453423 0.3318641  4.197750  0.043
tar_load(c(modfixdens12, modfdsel, modrdsel))
jtools::export_summs(
  modfixdens12,
  modfdsel,
  modrdsel,
  model.names = c("Density ~ Complexity", "Density ~ Complexity + Environment", "Density ~ Complexity + Environment + (1|Site)"),
  statistics = c(
    N = "nobs",
    DF = "df.residual",
    AIC = "AIC",
    BIC = "BIC",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
    coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Mean current velocity" = "Current_mean",
    "Depth",
    "Exposure" = "Fetch_max",
    "Gravel",
    "Organic Matter %" = "OM",
    "Mean Temperature" = "T_mean",
    "Year"
  ),
  scale = TRUE,
  error_pos = "right",
  bold_signif = .05,
  to.file = "docx",
  file.name = "outputs/summs_tot_dens.docx"
)
Density ~ ComplexityDensity ~ Complexity + EnvironmentDensity ~ Complexity + Environment + (1|Site)
Rhodolith complexity (PC1)0.19 ***(0.05)0.13 *  (0.05)0.03   (0.07)
Bed complexity (PC2)0.28 ***(0.05)-0.02    (0.06)0.11   (0.07)
PC1:PC2-0.35 ***(0.05)-0.19 ***(0.05)0.07   (0.06)
Mean current velocity           -0.27 ***(0.05)-0.51   (0.24)
Depth           -0.31 ***(0.06)0.09   (0.14)
Exposure           0.11    (0.07)0.15   (0.19)
Gravel           0.21 ***(0.05)0.05   (0.05)
Organic Matter %           0.18 ** (0.06)0.15 **(0.05)
Mean Temperature           0.09 *  (0.04)0.06   (0.03)
Year           0.04    (0.04)0.05   (0.03)
N348           348           348          
DF344.00        337.00        335.00       
AIC889.49        753.38        687.53       
BIC908.75        799.60        737.61       
Marginal R2                      0.24       
Conditional R20.26        0.52        0.70       
Adjusted R20.25        0.50                  
F39.46        36.03                  
p-value0.00        0.00                  
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.

By phylum

tar_load(c(modfdann, modrdann, modfdart, modrdart, modfdmol, modrdmol, modfdepi, modrdepi, modfdinf, modrdinf, modfdint, modrdint))

jtools::export_summs(
  modfdann,
  modrdann,
  modfdart,
  modrdart,
  modfdmol,
  modrdmol,
  model.names = c("Annelida LM","Annelida LMM", "Arthropoda LM", "Arthropoda LMM", "Mollusca LM", "Mollusca LMM"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Mean current velocity" = "Current_mean",
    "Depth",
    "Exposure" = "Fetch_max",
    "Gravel",
    "Organic Matter %" = "OM",
    "Mean Temperature" = "T_mean",
    "Year"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)
Annelida LMAnnelida LMMArthropoda LMArthropoda LMMMollusca LMMollusca LMM
Rhodolith complexity (PC1)-0.04    (0.05)0.05    (0.07)0.17 ** (0.06)0.07  (0.07)-0.26 ***(0.05)-0.14   (0.07)
Bed complexity (PC2)0.17 ** (0.06)0.14 *  (0.07)0.10    (0.06)0.04  (0.07)0.18 ***(0.05)0.11   (0.07)
PC1:PC20.01    (0.04)0.13 *  (0.06)-0.24 ***(0.05)0.08  (0.05)-0.13 ** (0.05)-0.01   (0.06)
Mean current velocity                      -0.28 ***(0.06)-0.09  (0.22)                     
Depth-0.29 ***(0.05)-0.21    (0.10)                    -0.46 ***(0.06)-0.32 **(0.11)
Exposure-0.25 ***(0.06)-0.38 ** (0.12)0.20 ** (0.07)0.20  (0.19)                     
Gravel                      0.28 ***(0.05)0.02  (0.04)0.15 ** (0.05)0.03   (0.05)
Organic Matter %0.09    (0.06)           0.22 ***(0.06)0.11 *(0.05)0.17 ***(0.05)0.15 **(0.06)
Mean Temperature0.09 *  (0.04)0.08 *  (0.04)                    0.09 *  (0.04)0.09 * (0.04)
Year0.14 ***(0.04)0.16 ***(0.03)0.06    (0.04)0.07 *(0.03)-0.06    (0.04)-0.06   (0.04)
N348           348           348           348         348           348          
Residual DF338.00        338.00        339.00        337.00      339.00        337.00       
Marginal R2           0.50                   0.03                 0.28       
Conditional R20.60        0.60        0.40        0.74      0.53        0.47       
Adjusted R20.59                   0.39                 0.52                  
F55.35                   28.65                 47.92                  
p-value0.00                   0.00                 0.00                  
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(anndenspc1)
tar_read(anndenspc2) 
tar_read(artdenspc1)
tar_read(artdenspc2)
tar_read(moldenspc1)
tar_read(moldenspc2)
Main phyla density as a function of HCMain phyla density as a function of HCMain phyla density as a function of HCMain phyla density as a function of HCMain phyla density as a function of HCMain phyla density as a function of HC

Main phyla density as a function of HC

By compartment

jtools::export_summs(
  modfdepi,
  modrdepi,
  modfdinf,
  modrdinf,
  modfdint,
  modrdint,
  model.names = c("Epifauna LM", "Epifauna LMM", "Infauna LM", "Infauna LMM", "Interstitial fauna LM", "Interstitial fauna LMM"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Mean current velocity" = "Current_mean",
    "Depth",
    "Exposure" = "Fetch_max",
    "Gravel",
    "Organic Matter %" = "OM",
    "Mud",
    "Mean Temperature" = "T_mean",
    "Temperature Variation (sd)" = "T_sd",
    "Year"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)
Epifauna LMEpifauna LMMInfauna LMInfauna LMMInterstitial fauna LMInterstitial fauna LMM
Rhodolith complexity (PC1)0.17 ** (0.06)0.07    (0.07)-0.00    (0.05)0.02   (0.08)-0.16 ** (0.05)-0.14 *  (0.07)
Bed complexity (PC2)0.13 *  (0.06)0.09    (0.07)0.23 ***(0.06)0.09   (0.08)0.21 ***(0.05)0.25 ***(0.07)
PC1:PC2-0.21 ***(0.05)0.05    (0.06)0.03    (0.05)0.16 **(0.06)-0.08    (0.04)0.01    (0.05)
Mean current velocity-0.27 ***(0.05)-0.24    (0.21)                                           
Depth                      -0.47 ***(0.06)-0.19   (0.12)                      
Exposure0.19 ** (0.07)0.15    (0.19)                     -0.32 ***(0.06)-0.34 ** (0.11)
Gravel0.34 ***(0.05)0.06    (0.05)-0.19 ***(0.05)-0.09   (0.05)0.28 ***(0.04)0.17 ** (0.05)
Organic Matter %0.27 ***(0.06)0.17 ***(0.05)                     0.07    (0.05)0.07    (0.05)
Mud                                           0.37 ***(0.06)0.25 ***(0.06)
Mean Temperature0.06    (0.04)0.04    (0.03)0.09    (0.04)0.06   (0.04)                      
Temperature Variation (sd)                      -0.05    (0.05)-0.01   (0.06)                      
Year0.00    (0.04)0.02    (0.03)0.09 *  (0.04)0.10 * (0.04)0.14 ***(0.04)0.14 ***(0.03)
N348           348           348           348          348           348           
Residual DF338.00        336.00        339.00        337.00       339.00        337.00        
Marginal R2           0.10                   0.13                  0.53        
Conditional R20.45        0.66        0.43        0.52       0.60        0.63        
F30.68                   31.92                  63.38                   
p-value0.00                   0.00                  0.00                   
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(epidenspc1)
tar_read(epidenspc2) 
tar_read(infdenspc1)
tar_read(infdenspc2)
tar_read(intdenspc1)
tar_read(intdenspc2)
Density of macrofauna as a function of HC at each compartmentDensity of macrofauna as a function of HC at each compartmentDensity of macrofauna as a function of HC at each compartmentDensity of macrofauna as a function of HC at each compartmentDensity of macrofauna as a function of HC at each compartmentDensity of macrofauna as a function of HC at each compartment

Density of macrofauna as a function of HC at each compartment 

tar_read(psdens)
tar_read(psdtrait)
Figure 2 (main text)Figure 2 (main text)

Figure 2 (main text)

Observed richness

Check if the relationships between HC and species richness are the same for the total macrofauna as well as for the main phyla and the main compartments (different positions in sediment) separately.

Total macrofauna

Now we first forward select environmental variables and then add sites as a random factor to control for pseudo-replication and other site-related variables that were not taken into account.

tar_read(totrichpc1)
tar_read(figure3)
Figure 3. OLS linear regression of Total Macrofauna Richness as a function of HCFigure 3. OLS linear regression of Total Macrofauna Richness as a function of HC

Figure 3. OLS linear regression of Total Macrofauna Richness as a function of HC 

tar_read(richselR2)
##   variables order          R2     R2Cum  AdjR2Cum          F pvalue
## 1     Depth     7 0.329480820 0.3294808 0.3275429 170.018051  0.001
## 2 Fetch_max     8 0.049649871 0.3791307 0.3755314  27.589068  0.001
## 3    T_mean     4 0.019017288 0.3981480 0.3928993  10.869694  0.002
## 4 PC2_score    10 0.009087146 0.4072351 0.4003224   5.258225  0.017
## 5        OM     3 0.013428474 0.4206636 0.4121938   7.927239  0.005
tar_load(c(modfixrich12, modfrsel, modrrsel))

jtools::export_summs(
  modfixrich12,
  modfrsel,
  modrrsel,
  model.names = c("S ~ Complexity", "S ~ Complexity + Environment", "S ~ Complexity + Environment + (1|Site)"),
  statistics = c(
    N = "nobs",
    DF = "df.residual",
    AIC = "AIC",
    BIC = "BIC",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
  coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Depth",
    "Exposure" = "Fetch_max",
    "Mean bottom temperature" = "T_mean",
    "Organic Matter %" = "OM",
    "Year"
  ),
  scale = TRUE,
  error_pos = "right",
  to.file = "docx",
  bold_signif = .05,
  file.name = "outputs/summs_tot_rich.docx"
)
S ~ ComplexityS ~ Complexity + EnvironmentS ~ Complexity + Environment + (1|Site)
Rhodolith complexity (PC1)2.42 *  (0.97)-0.47    (1.05)1.71   (1.60)
Bed complexity (PC2)10.74 ***(0.93)4.41 ** (1.36)3.06   (1.70)
PC1:PC2-1.01    (1.00)-1.08    (0.96)1.96   (1.34)
Depth           -5.56 ***(1.27)-2.77   (2.58)
Exposure           -2.97 *  (1.42)-1.73   (3.18)
Mean bottom temperature           1.97 *  (0.90)1.39   (0.87)
Organic Matter %           3.09 *  (1.25)3.98 **(1.25)
Year           1.68    (0.88)1.76 * (0.84)
N348           348           348          
DF344.00        339.00        337.00       
AIC2973.86        2909.46        2877.62       
BIC2993.12        2947.98        2920.00       
Marginal R2                      0.27       
Conditional R20.29        0.43        0.44       
Adjusted R20.29        0.42                  
F47.43        31.81                  
p-value0.00        0.00                  
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.

You might notice that PC1 was not selected as it has a very small effect size and doesn’t contribute enough to the adjusted R2. I kept it in the final model anyways for comparison with the complexity-only model.

tar_read(psrich)
Figure 4 (main text)

Figure 4 (main text)

By phylum

tar_load(c(modrrann, modrrart, modrrmol, modrrepi, modrrinf, modrrint))

jtools::export_summs(
  modrrann,
  modrrart,
  modrrmol,
  model.names = c("Annelida", "Arthropoda", "Mollusca"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
  coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score"
  ),
  scale = TRUE,
  error_pos = "right",
  bold_signif = .05
)
AnnelidaArthropodaMollusca
Rhodolith complexity (PC1)2.16 **(0.79)0.62(0.54)-0.77(0.58)
Bed complexity (PC2)1.87 * (0.81)0.58(0.56)0.25(0.60)
PC1:PC22.02 **(0.67)0.42(0.44)-0.75(0.48)
N348          348       348       
Residual DF342.00       342.00    342.00    
Marginal R20.12       0.02    0.02    
Conditional R20.33       0.52    0.48    
p-value                        
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(annrichpc1)
tar_read(annrichpc2) 
tar_read(artrichpc1)
tar_read(artrichpc2)
tar_read(molrichpc1)
tar_read(molrichpc2)
Main phyla richness as a function of HCMain phyla richness as a function of HCMain phyla richness as a function of HCMain phyla richness as a function of HCMain phyla richness as a function of HCMain phyla richness as a function of HC

Main phyla richness as a function of HC

By compartment

jtools::export_summs(
  modrrepi,
  modrrinf,
  modrrint,
  model.names = c("Epifauna", "Infauna", "Interstitial fauna"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
    coefs = c(
    "Rhodolith complexity (PC1)" = "PC1_score", 
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "PC1_score:PC2_score"
  ),
  scale = TRUE,
  error_pos = "right",
  bold_signif = .05
)
EpifaunaInfaunaInterstitial fauna
Rhodolith complexity (PC1)0.84(0.92)1.26  (0.69)0.31   (0.36)
Bed complexity (PC2)1.09(0.96)0.21  (0.71)1.13 **(0.37)
PC1:PC20.67(0.76)1.44 *(0.57)0.17   (0.31)
N348       348         347          
Residual DF342.00    342.00      341.00       
Marginal R20.02    0.03      0.08       
Conditional R20.51    0.48      0.30       
p-value                          
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(epirichpc1)
tar_read(epirichpc2) 
tar_read(infrichpc1)
tar_read(infrichpc2)
tar_read(intrichpc1)
tar_read(intrichpc2)
Species richness as a function of HC at each faunal compartmentSpecies richness as a function of HC at each faunal compartmentSpecies richness as a function of HC at each faunal compartmentSpecies richness as a function of HC at each faunal compartmentSpecies richness as a function of HC at each faunal compartmentSpecies richness as a function of HC at each faunal compartment

Species richness as a function of HC at each faunal compartment

As we can see, the marginal R2 for the models containing only a subset of the macrofauna are quite low. However, although the effects of HC are less evident in these cases, The pattern for most subsets is very similar to the general patter, with very slight positive effect of rhodolith complexity and a more important positive effect of bed complexity. Nevertheless Rhodolith complexity seems to play a more important role in driving annelid richness than the other main phyla.

N2

tar_read(psN2)
Regression coefficients of LMs and LMM of Simpson's inverse diversity (N2) as a function of HC and physical environmental variables.

Regression coefficients of LMs and LMM of Simpson’s inverse diversity (N2) as a function of HC and physical environmental variables.

J

tar_read(psJ)
Regression coefficients of LMs and LMM of Pielou's evenness (J) as a function of HC and physical

Regression coefficients of LMs and LMM of Pielou’s evenness (J) as a function of HC and physical

Regional patterns of biodiversity

PCA

tar_read(spfauna)
Screeplot of the PCA of log-chord transformed density data of total macrofauna

Screeplot of the PCA of log-chord transformed density data of total macrofauna

Final PCA showing only the species with a godness of fit \(\geq\) .3

tar_read(faunapca)
Figure 5 (main text)

Figure 5 (main text)

Cluster

tar_read(faunaclust)
Ward's cluster of median complexity values at the site level

Ward’s cluster of median complexity values at the site level

Effects of HC on regional diversity patterns

Redundancy analysis

After variable selection, we add sites as a factor to understand if we describe them well with the environmental variables we chose. We can look at the first two axes to have an idea of the relationships between the different explanatory variables.

tar_read(rdafsel)
##            variables order          R2     R2Cum   AdjR2Cum         F pvalue
## 1                Mud     2 0.087006599 0.0870066 0.08436789 32.973166  0.001
## 2              Depth     8 0.052115736 0.1391223 0.13413174 20.885579  0.001
## 3       Current_mean     7 0.052409189 0.1915315 0.18448093 22.299894  0.001
## 4  Branching_density    10 0.037531459 0.2290630 0.22007246 16.698239  0.001
## 5              D_bin    11 0.027960583 0.2570236 0.24616134 12.870555  0.001
## 6          Fetch_max     9 0.021665334 0.2786889 0.26599721 10.242292  0.001
## 7               T_sd     6 0.017883207 0.2965721 0.28208977  8.643800  0.001
## 8         Sphericity    15 0.013192976 0.3097651 0.29347635  6.479560  0.001
## 9               Year     1 0.011973556 0.3217386 0.30367842  5.966818  0.001
## 10               DR3    16 0.010094019 0.3318327 0.31200573  5.091067  0.001
## 11            Gravel     3 0.009114600 0.3409473 0.31937113  4.646829  0.001
## 12            T_mean     5 0.008733095 0.3496804 0.32638532  4.498691  0.001
## 13            D_gray    12 0.007959468 0.3576398 0.33263778  4.138585  0.001
## 14                 L    13 0.005985208 0.3636250 0.33687053  3.131918  0.001
## 15     Total_Density    14 0.005467386 0.3690924 0.34058755  2.877081  0.001
## 16                OM     4 0.004910747 0.3740032 0.34374350  2.596590  0.001
tar_load(rdasite)
tar_read(triplotrda)
RDA triplot showing the variability in species composition related to physical environmental variables and habitat complexity. The first two axes are shown.

RDA triplot showing the variability in species composition related to physical environmental variables and habitat complexity. The first two axes are shown. 

RsquareAdj(rdasite)$adj.r.squared
## [1] 0.4297523

We check the model’s validity. Only 99 permutations were chosen for the validation by axis and by term here as they are really time consuming but the results are the same with 999 permutations. (Output is collapsed, click on code to see)

tar_read(aovsite)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = bcdens ~ Mud + Depth + Current_mean + Branching_density + D_bin + Fetch_max + T_sd + Sphericity + T_mean + DR3 + Gravel + D_gray + Total_Density + L + OM + Year + Site, data = envcomp)
##           Df Variance     F Pr(>F)    
## Model     25  0.28809 11.46  0.001 ***
## Residual 322  0.32377                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tar_read(aovaxis)
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 99
## 
## Model: rda(formula = bcdens ~ Mud + Depth + Current_mean + Branching_density + D_bin + Fetch_max + T_sd + Sphericity + T_mean + DR3 + Gravel + D_gray + Total_Density + L + OM + Year + Site, data = envcomp)
##           Df Variance       F Pr(>F)   
## RDA1       1  0.07622 75.8006   0.01 **
## RDA2       1  0.04629 46.0328   0.01 **
## RDA3       1  0.04255 42.3172   0.01 **
## RDA4       1  0.02862 28.4608   0.01 **
## RDA5       1  0.01877 18.6709   0.01 **
## RDA6       1  0.01353 13.4585   0.01 **
## RDA7       1  0.01227 12.2028   0.01 **
## RDA8       1  0.00965  9.5928   0.01 **
## RDA9       1  0.00784  7.8011   0.01 **
## RDA10      1  0.00720  7.1584   0.01 **
## RDA11      1  0.00589  5.8598   0.01 **
## RDA12      1  0.00308  3.0642   0.01 **
## RDA13      1  0.00267  2.6539   0.01 **
## RDA14      1  0.00212  2.1037   0.17   
## RDA15      1  0.00180  1.7924   0.51   
## RDA16      1  0.00155  1.5436   0.86   
## RDA17      1  0.00133  1.3219   0.99   
## RDA18      1  0.00129  1.2831   0.99   
## RDA19      1  0.00106  1.0509   1.00   
## RDA20      1  0.00091  0.9030   1.00   
## RDA21      1  0.00084  0.8355   1.00   
## RDA22      1  0.00082  0.8172   1.00   
## RDA23      1  0.00073  0.7300   1.00   
## RDA24      1  0.00059  0.5858   1.00   
## RDA25      1  0.00047  0.4668   1.00   
## Residual 322  0.32377                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tar_read(aovterm)
## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 99
## 
## Model: rda(formula = bcdens ~ Mud + Depth + Current_mean + Branching_density + D_bin + Fetch_max + T_sd + Sphericity + T_mean + DR3 + Gravel + D_gray + Total_Density + L + OM + Year + Site, data = envcomp)
##                    Df Variance       F Pr(>F)   
## Mud                 1  0.05324 52.9442   0.01 **
## Depth               1  0.03189 31.7128   0.01 **
## Current_mean        1  0.03207 31.8914   0.01 **
## Branching_density   1  0.02296 22.8382   0.01 **
## D_bin               1  0.01711 17.0142   0.01 **
## Fetch_max           1  0.01326 13.1835   0.01 **
## T_sd                1  0.01094 10.8821   0.01 **
## Sphericity          1  0.00807  8.0280   0.01 **
## T_mean              1  0.00695  6.9100   0.01 **
## DR3                 1  0.00610  6.0709   0.01 **
## Gravel              1  0.00556  5.5328   0.01 **
## D_gray              1  0.00487  4.8401   0.01 **
## Total_Density       1  0.00334  3.3223   0.01 **
## L                   1  0.00366  3.6403   0.01 **
## OM                  1  0.00318  3.1599   0.01 **
## Year                1  0.00564  5.6130   0.01 **
## Site                9  0.05925  6.5471   0.01 **
## Residual          322  0.32377                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hierarchical partitioning

tar_read(upsetmain)
Figure 6 (main text)

Figure 6 (main text) 

tar_read(upsetsep)
Hierarchical partitioning separating Rhodolith and Bed complexity

Hierarchical partitioning separating Rhodolith and Bed complexity

Community stability

tar_read(bdplot)
Figure 7 (main text)

Figure 7 (main text) 

tar_read(replplot)
Relative importance of the Replacement component of temporal BD

Relative importance of the Replacement component of temporal BD

Effects of HC on community stability

Total Macrofauna

tar_load(c(modbdtot, modrdftot, modrptot))
jtools::export_summs(
  modbdtot,
  modrptot,
  modrdftot,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)
BDtotalReplacementRichness Difference
Rhodolith complexity (PC1)-0.34 * (0.15)-0.29  (0.18)-0.18(0.19)
Rhodolith complexity (PC1)^2-0.52 **(0.16)-0.38  (0.19)-0.36(0.21)
Bed complexity (PC2)-0.51 **(0.15)-0.46 *(0.17)-0.24(0.19)
PC1:PC20.15   (0.18)0.15  (0.21)0.05(0.22)
PC1ˆ2:PC2-0.23   (0.21)-0.30  (0.24)0.01(0.27)
N30          30         30       
Residual DF24.00       24.00      24.00    
R20.49       0.33      0.20    
Adjusted R20.38       0.20      0.03    
F4.58       2.41      1.21    
p-value0.00       0.07      0.34    
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(bdtotpc1)
Effects of HC on community stability

Effects of HC on community stability 

tar_read(bdtotpc2)
Effects of HC on community stability

Effects of HC on community stability 

tar_read(rdtotpc1)
Effects of HC on community stability

Effects of HC on community stability 

tar_read(rdtotpc2)
Effects of HC on community stability

Effects of HC on community stability 

tar_read(rptotpc1)
Effects of HC on community stability

Effects of HC on community stability 

tar_read(rptotpc2)
Effects of HC on community stability

Effects of HC on community stability

By phylum

Annelida

tar_load(c(modbdann, modrdfann, modrpann))
jtools::export_summs(
  modbdann,
  modrpann,
  modrdfann,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)
BDtotalReplacementRichness Difference
Rhodolith complexity (PC1)-0.51 ** (0.14)-0.60 ***(0.14)-0.04(0.19)
Rhodolith complexity (PC1)^2-0.51 ** (0.15)-0.33 *  (0.15)-0.40(0.20)
Bed complexity (PC2)-0.52 ***(0.14)-0.44 ** (0.14)-0.28(0.19)
PC1:PC20.05    (0.16)0.11    (0.17)-0.06(0.22)
PC1ˆ2:PC2-0.24    (0.19)-0.42 *  (0.20)0.16(0.26)
N30           30           30       
Residual DF24.00        24.00        24.00    
R20.59        0.56        0.24    
Adjusted R20.51        0.47        0.08    
F6.93        6.06        1.48    
p-value0.00        0.00        0.23    
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(bdannpc1)
Effects of HC on annelid community stability

Effects of HC on annelid community stability 

tar_read(bdannpc2)
Effects of HC on annelid community stability

Effects of HC on annelid community stability 

tar_read(rdannpc1)
Effects of HC on annelid community stability

Effects of HC on annelid community stability 

tar_read(rdannpc2)
Effects of HC on annelid community stability

Effects of HC on annelid community stability 

tar_read(rpannpc1)
Effects of HC on annelid community stability

Effects of HC on annelid community stability 

tar_read(rpannpc2)
Effects of HC on annelid community stability

Effects of HC on annelid community stability

Arthropoda

tar_load(c(modbdart, modrdfart, modrpart))
jtools::export_summs(
  modbdart,
  modrpart,
  modrdfart,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)
BDtotalReplacementRichness Difference
Rhodolith complexity (PC1)-0.24  (0.17)-0.20  (0.18)-0.18(0.17)
Rhodolith complexity (PC1)^2-0.39 *(0.19)-0.34  (0.20)-0.29(0.19)
Bed complexity (PC2)-0.47 *(0.17)-0.46 *(0.18)-0.31(0.17)
PC1:PC20.21  (0.20)-0.14  (0.21)0.41(0.20)
PC1ˆ2:PC2-0.19  (0.24)-0.20  (0.25)-0.12(0.24)
N30         30         30       
Residual DF24.00      24.00      24.00    
R20.35      0.28      0.35    
Adjusted R20.22      0.13      0.21    
F2.59      1.84      2.53    
p-value0.05      0.14      0.06    
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(bdartpc1)
Effects of HC on arthropod community stability

Effects of HC on arthropod community stability 

tar_read(bdartpc2)
Effects of HC on arthropod community stability

Effects of HC on arthropod community stability 

tar_read(rdartpc1)
Effects of HC on arthropod community stability

Effects of HC on arthropod community stability 

tar_read(rdartpc2)
Effects of HC on arthropod community stability

Effects of HC on arthropod community stability 

tar_read(rpartpc1)
Effects of HC on arthropod community stability

Effects of HC on arthropod community stability 

tar_read(rpartpc2)
Effects of HC on arthropod community stability

Effects of HC on arthropod community stability

Mollusca

tar_load(c(modbdmol, modrdfmol, modrpmol))
jtools::export_summs(
  modbdmol,
  modrpmol,
  modrdfmol,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)
BDtotalReplacementRichness Difference
Rhodolith complexity (PC1)0.03   (0.17)0.11(0.20)-0.11  (0.19)
Rhodolith complexity (PC1)^2-0.50 * (0.18)-0.33(0.21)-0.25  (0.20)
Bed complexity (PC2)-0.52 **(0.17)-0.25(0.19)-0.40 *(0.19)
PC1:PC20.10   (0.20)0.18(0.23)-0.13  (0.22)
PC1ˆ2:PC2-0.20   (0.24)-0.03(0.27)-0.26  (0.27)
N30          30       30         
Residual DF24.00       24.00    24.00      
R20.36       0.17    0.21      
Adjusted R20.23       -0.01    0.05      
F2.69       0.96    1.28      
p-value0.05       0.46    0.31      
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(bdmolpc1)
Effects of HC on mollusc community stability

Effects of HC on mollusc community stability 

tar_read(bdmolpc2)
Effects of HC on mollusc community stability

Effects of HC on mollusc community stability 

tar_read(rdmolpc1)
Effects of HC on mollusc community stability

Effects of HC on mollusc community stability 

tar_read(rdmolpc2)
Effects of HC on mollusc community stability

Effects of HC on mollusc community stability 

tar_read(rpmolpc1)
Effects of HC on mollusc community stability

Effects of HC on mollusc community stability 

tar_read(rpmolpc2)
Effects of HC on mollusc community stability

Effects of HC on mollusc community stability

By compartment

Epifauna

tar_load(c(modbdepi, modrdfepi, modrpepi))
jtools::export_summs(
  modbdepi,
  modrpepi,
  modrdfepi,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)
BDtotalReplacementRichness Difference
Rhodolith complexity (PC1)-0.32  (0.17)-0.19(0.20)-0.28(0.18)
Rhodolith complexity (PC1)^2-0.45 *(0.18)-0.30(0.21)-0.38(0.20)
Bed complexity (PC2)-0.33  (0.17)-0.23(0.20)-0.26(0.18)
PC1:PC20.23  (0.20)0.16(0.23)0.18(0.22)
PC1ˆ2:PC2-0.20  (0.24)-0.15(0.27)-0.15(0.26)
N30         30       30       
Residual DF24.00      24.00    24.00    
R20.37      0.16    0.26    
Adjusted R20.24      -0.02    0.11    
F2.84      0.89    1.72    
p-value0.04      0.51    0.17    
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(bdepipc1)
Effects of HC on epifaunal community stability

Effects of HC on epifaunal community stability 

tar_read(bdepipc2)
Effects of HC on epifaunal community stability

Effects of HC on epifaunal community stability 

tar_read(rdepipc1)
Effects of HC on epifaunal community stability

Effects of HC on epifaunal community stability 

tar_read(rdepipc2)
Effects of HC on epifaunal community stability

Effects of HC on epifaunal community stability 

tar_read(rpepipc1)
Effects of HC on epifaunal community stability

Effects of HC on epifaunal community stability 

tar_read(rpepipc2)
Effects of HC on epifaunal community stability

Effects of HC on epifaunal community stability

Infauna

tar_load(c(modbdinf, modrdfinf, modrpinf))
jtools::export_summs(
  modbdinf,
  modrpinf,
  modrdfinf,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)
BDtotalReplacementRichness Difference
Rhodolith complexity (PC1)-0.43 ** (0.13)-0.38 * (0.16)-0.14(0.19)
Rhodolith complexity (PC1)^2-0.56 ***(0.14)-0.49 **(0.17)-0.18(0.21)
Bed complexity (PC2)-0.65 ***(0.13)-0.49 **(0.16)-0.36(0.19)
PC1:PC2-0.09    (0.15)-0.00   (0.19)-0.16(0.23)
PC1ˆ2:PC2-0.28    (0.18)-0.41   (0.22)0.19(0.27)
N30           30          30       
Residual DF24.00        24.00       24.00    
R20.65        0.43       0.19    
Adjusted R20.57        0.31       0.02    
F8.82        3.66       1.14    
p-value0.00        0.01       0.36    
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(bdinfpc1)
Effects of HC on infaunal community stability

Effects of HC on infaunal community stability 

tar_read(bdinfpc2)
Effects of HC on infaunal community stability

Effects of HC on infaunal community stability 

tar_read(rdinfpc1)
Effects of HC on infaunal community stability

Effects of HC on infaunal community stability 

tar_read(rdinfpc2)
Effects of HC on infaunal community stability

Effects of HC on infaunal community stability 

tar_read(rpinfpc1)
Effects of HC on infaunal community stability

Effects of HC on infaunal community stability 

tar_read(rpinfpc2)
Effects of HC on infaunal community stability

Effects of HC on infaunal community stability

Interstitial fauna

tar_load(c(modbdint, modrdfint, modrpint))
jtools::export_summs(
  modbdint,
  modrpint,
  modrdfint,
  model.names = c("BDtotal", "Replacement", "Richness Difference"),
  statistics = c(
    N = "nobs",
    `Residual DF` = "df.residual",
    `R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
   coefs = c(
    "Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
    "Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
    "Bed complexity (PC2)" = "PC2_score",
    "PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
    "PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
  ),
  scale = TRUE,
  transform.response = TRUE,
  error_pos = "right",
  bold_signif = .05
)
BDtotalReplacementRichness Difference
Rhodolith complexity (PC1)0.03   (0.17)-0.20  (0.18)0.24(0.18)
Rhodolith complexity (PC1)^2-0.39 * (0.18)-0.01  (0.19)-0.40(0.20)
Bed complexity (PC2)-0.56 **(0.17)-0.45 *(0.17)-0.13(0.18)
PC1:PC20.27   (0.19)0.36  (0.21)-0.09(0.21)
PC1ˆ2:PC2-0.08   (0.23)-0.29  (0.24)0.21(0.25)
N30          30         30       
Residual DF24.00       24.00      24.00    
R20.40       0.34      0.28    
Adjusted R20.27       0.20      0.13    
F3.19       2.42      1.85    
p-value0.02       0.06      0.14    
All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05.
tar_read(bdintpc1)
Effects of HC on intersitial community stability

Effects of HC on intersitial community stability 

tar_read(bdintpc2)
Effects of HC on intersitial community stability

Effects of HC on intersitial community stability 

tar_read(rdintpc1)
Effects of HC on intersitial community stability

Effects of HC on intersitial community stability 

tar_read(rdintpc2)
Effects of HC on intersitial community stability

Effects of HC on intersitial community stability 

tar_read(rpintpc1)
Effects of HC on intersitial community stability

Effects of HC on intersitial community stability 

tar_read(rpintpc2)
Effects of HC on intersitial community stability

Effects of HC on intersitial community stability 

tar_read(bdtrait)
tar_read(bdphyl)
Figure 8 (main text)Figure 8 (main text)

Figure 8 (main text)

Box, G. E. P., and D. R. Cox. 1964. “An Analysis of Transformations.” Journal of the Royal Statistical Society. Series B (Methodological) 26 (2): 211–52. https://www.jstor.org/stable/2984418.
Casajus N. 2022. Rcompendium: An r Package to Create a Package or Research Compendium Structure. https://github.com/FRBCesab/rcompendium.
Landau, William Michael. 2021. “The Targets r Package: A Dynamic Make-Like Function-Oriented Pipeline Toolkit for Reproducibility and High-Performance Computing.” Journal of Open Source Software 6 (57): 2959. https://doi.org/10.21105/joss.02959.